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^ ABSTRACT 

This paper presents numerical solutions of particular Tolman models approx- 
£NJ imating a fractal behaviour along the past light cone. The initial conditions of 

the numerical problem are discussed and the algorithm used to carry out the nu- 
j — merical integrations is presented. It was found that the numerical solutions are 

stiff across the flat-curved interface necessary to obtain the initial conditions of 
the problem. The spatially homogeneous Friedmann models are treated as special 
m ^ cases of the Tolman solution and solved numerically. Extending the results of pa- 

per II on the Einstein-de Sitter model, to the K = ±1 models, it was found that 

H 

the open and closed Friedmann models also do not appear to remain homogeneous 
along the backward null cone, with a vanishing volume (average) density as one 
approaches the big bang singularity hypersurface. Fractal solutions, that is, solu- 
tions representing an averaged and smoothed-out single fractal, were obtained in 



1 Present address: Departamento de Astroffsica, Observatorio Nacional - CNPq, Rua General Jose 
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all three classes of the Tolman metric, but only the hyperbolic ones were found 
to be in agreement with observations, meaning that a possible Friedmann back- 
ground universe would have to be an open one. The best fractal metric obtained 
through numerical simulations is also analysed in terms of evolution, homothetic 
self-similarity, comparison with the respective spatially homogeneous case and the 
fitting problem in cosmology. The paper finishes with a discussion on some objec- 
tions raised by some authors against a fractal cosmology. 



1 Introduction 

This is the third paper of a series that studies a relativistic cosmology modelling the 
relativistic generalization of the single fractal Newtonian model advanced by Pietronero 
(1987; see also Coleman & Pietronero 1992), in which the galactic clustering problem is 
studied by assuming that the large-scale structure of the universe can be described as 
being a self-similar fractal system]^] 

In Ribeiro (1992a, hereafter paper I) I argued that the recent all sky redshift sur- 
veys (de Lapparent, Geller & Huchra 1986; Saunders et al. 1991) present observations 
consistent both with the old Charlier hypothesis of hierarchical clustering and with frac- 
tals, where the latter is in essence a more precise conceptualization of the scaling idea 
implicit in the hierarchical clustering hypothesis. In paper I Pietronero's (1987) basic 
hypotheses were assumed in order to propose similar ones in a relativistic context, and I 
obtained observational relations compatible with fractals in Tolman's spacetime and de- 
vised a numerical strategy for finding particular Tolman solutions representing a fractal 
behaviour along the backward null cone. In Ribeiro (1992b, hereafter paper II) I studied 
analytically the Einstein-de Sitter model in the context of the theory developed in pa- 
per I. By treating the Einstein-de Sitter model as a special case of Tolman's spacetime, 
I found that it does not appear to remain homogeneous along the past null geodesic, 
has a volume (average) density which vanishes asymptotically and that it also shows 

2 The models investigated in this series of papers are in the realm of classical cosmology. No 
hypotheses concerning inflationary cosmology have been considered. 
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no single fractal features along the backward null cone. The apparent inhomogeneity of 
the Einstein-de Sitter model is explained by the fact that densities measured along the 
geodesic go through different hypersurfaces of constant t, where each one has a different 
value for the proper density. 

This paper continues the study of these cosmologies and presents relativistic fractal 
solutions obtained by following the numerical simulation strategy already devised in 
paper I. By fractal solutions I mean solutions where the fractal system is smoothed-out 
and its average density follows the de Vaucouleurs density power-law. These solutions 
represent fractal behaviour along the backward null cone and they were obtained for 
all three types of Tolman dust models, namely, elliptic, parabolic and hyperbolic. By 
analysing these solutions we conclude that the only ones with features that may represent 
real astronomical observations are of hyperbolic type. As we are studying the Tolman 
spacetime as a region (or regions) possibly surrounded by a Friedmann universe (see 
paper I), if we adopt the fitting condition approach for interpreting the match between 
the two spacetimes (Ellis & Stoeger 1987; Ellis & Jaklitsch 1989) we then conclude that 
the Friedmann background required is also hyperbolic and we would be living in an open, 
ever expanding, universe. 

This paper is organized as follows. In §2 is shown a summary of the observational 
relations developed in papers I and II plus some minor extensions which will be necessary 
here, and in §3 I discuss the initial conditions and the algorithm used to find Tolman 
numerical solutions along the past light cone. In §4 I present the numerical solutions 
for the spatially homogeneous special cases and §5 shows fractal solutions for all types 
of Tolman models. §6 discusses these fractal solutions in terms of comparison with the 
spatially homogeneous cases, fitting condition, evolution of the most realistic fractal 
model in terms of real observations and relations with homothetic self-similarity. In this 
section I also express my criticism of criticisms of fractal cosmology. I finish in §7 with 
the conclusions on this and the preceding papers. 
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2 Observational Relations in Tolman's Dust Models 



The Tolman (1934) metric for the motion of spherically symmetric dust is 

dS 2 = dt 2 - ^dr 2 - R 2 (d6 2 + sin 2 6d<p 2 ). (1) 

The Einstein field equations (with c = G — 1 and A = 0) for metric ([!]) reduce to a 
single equation 

2RR 2 + 2R(1 - f 2 ) = F, (2) 

where a dot means d/dt and a prime means d/dr; f(r) and F(r) are two arbitrary 
functions and R(r, t) > 0. The analytical solutions of equation ^ are divided into three 
categories according to the values of the function f(r). For f 2 < 1 we have the elliptic 
class and equation ^ has the following solution: 

F(l-cos29) 

R= 4|/2-l | ' (3) 

where is given by 

F(2e-sin26) 
4| f 2 - 1 

For / 2 = 1 we have the parabolic class and the solution of equation ^ may be written 

as 

R= I(9F) 1/3 (t + /3) 2/3 . (5) 

Finally the hyperbolic class is given when f 2 > 1, and equation ^ is solved by 

_ F( C03 h2e-l) 

4(/»-l) ' (6) 

where 

F(sinh29-2e) 

+ 4(/2-l) 3 / 2 • (7j 

Note that one Tolman solution can have regions of differing classes, depending on the 
values of the function f(r). 
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The derivatives R', R, R' for the three cases above form long algebraic expressions 
and I shall not show them here. They can be found in paper I. The function (3{r) is a 
third arbitrary function which gives the local time passed since the singularity surface, 
that is, since the big bang. In the elliptic and hyperbolic models, it is necessary to solve 
equations Q and ^ in terms of 9 so that specific values of R and its derivatives are 
evaluated. In order to do so numerically, we need to find the interval where the root lies. 
It is not difficult to show (see paper I) that in the elliptic models 

i(t + p)\ f - l | 3/2 - l < 26 < i(f + (3)\ f - 1 | 3/2 + 1, (8) 



and in the hyperbolic case 

< < 



|(t + /3)(/ 2 -l) 3/2 ^ /3 



(9) 



The local density in Tolman's spacetime is given by 



^=2im (10) 

and if we adopt the radius coordinate r as the parameter along the backward null cone 
(r > 0), we can then write the past radial null geodesic of metric ([I]) as 

The Friedmann metric is a special case of the metric ([I]) above, obtained when R = 
a(t) g(r), f = g' and (3 = constant. The usual Friedmann universe requires further that 
g = sinr, r, sinhr and F = b\ sin 3 r, |r 3 , 6 2 sinh 3 r, which are respectively the cases for 
K = +1, 0, —1. The positive constants b\ and hi are scaling factors necessary to make 
the density parameter Q equal to any value different from one in the open and closed 
models. It is also possible to show that F(r) gives the gravitational mass inside the 
comoving radius r and f(r) gives the total energy of the system also within r (see §6.1 
and paper II for further details). 

It was shown in paper II how the Hubble constant relates to the function j3(r) in a 
Einstein-de Sitter universe. It is important here to extend that result to the K = ±1 
Friedmann cases. By definition, H = R/R = a /a for constant j3(r). Considering the 
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especializations above, that permit us to get the Friedmann metric from the Tolman 
solution, and assuming t = as our "now", that is, t — being the time coordinate 
label for the present epoch, it is straightforward to show that the present value H for 
the Hubble constant in the closed Friedmann model is given by 

H = 4Sm29 ° 2 , (12) 
6 1 (l-cos26 ) 2 

where 0o is the solution of 

4/5 o = 6i(20 o -sin2e o ), (13) 

and /3(r) = (5 is a constant that gives the age of the universe. In the open Friedmann 
case we will have then 

Ho = 45inh29 ° a , (14) 
6 2 (cosh26o-l) 

and 

4/5o = 6 2 (sinh20o-2eo). (15) 

It will be necessary in further calculations to obtain the value of the cosmological 
density parameter Q in the Friedmann model at the present constant time hypersurface. 
By definition VL = p/ p c and p c = l/(Girf3o 2 ) at t — 0. Therefore, in the closed Friedmann 
model we have that 

72/V 

n ° = 7^i\ ( 16) 

(oi) (1 — cos2B ) 

where 0o is given by equation (13). In the open Friedmann model we have 

n ° = (6 2 ) 2 (cosh2e -l) 3 ' 



with O being the solution of equation ([15J). 

It was shown in paper I that in a Tolman universe the redshift may be written as 

l + z = (l-iy\ (18) 

where the function I(r) is the solution of the differential equation 
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The luminosity distance di and the cumulative number count N c are given by 

d, = R(l + zf, (20) 

where Mq is the average galactic rest mass (~ iO n M ). The volume V of the sphere 
which contains the sources and the volume density (average density) p v have the form 

4 



r = g7r(d,) a , 



Pv 



N C M G 
V ' 



(22) 



(23) 



The proposed relativistic version of Pietronero's (1987) generalized mass-length relation 
is given by 



N c = a{di) D (24) 
where a is a constant related to the lower cutoff of the fractal system and D is its fractal 



dimension. By substituting equations (22) and (24) into equation (23) we get the de 
Vaucouleurs' density power law 

3a Mq 



Pv 



4tt 
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D. 



(25) 



3 Numerical Analysis 



3.1 Initial Conditions 



The numerical task consists of finding solutions of the geodesic equation (11) and the 



equation (19) that enables the calculation of the redshift. Let us call the solution of the 



former t(r) and the solution of the latter I(r). As we are seeking results along the past 
light cone, at each step of the independent variable r iy we need to find ti first and then 
use it to find U {i = 1, 2, . . . , n). 

In paper I it was initially proposed to start the integration at t(0) = and 1(0) = 0, 
but realized that in this approach there would be problems in the elliptic and hyperbolic 
models because we would get [j type indeterminations in the function R(r, t) at r = 0. 
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In order to overcome this difficulty I then proposed a flat metric from r = till r = e (e 
small), beyond which the spacetime changed to a curved one. This hypothesis effectively 
replaces the initial conditions above by t(e) = —e, 1(e) = 0. As we shall see, these 
new initial conditions produce numerical results which are perfectly acceptable [j and 
numerically stable, but in this approach the formulation and solution of the interface 
between flat and curved spacetimes is left incomplete inasmuch as a full determination of 
this interface demands the solution of the junction conditions between the two regions. 
If this matching is not fully solved, discontinuities might arise across the joining surface, 
mainly in the function R(r, t) which determines the luminosity distance. As we are 
assuming e small (at least 10 -3 ) these possible discontinuities do not necessarily represent 
a hazard to the numerical solutions, but we need to have a way to check them. 

The match between flat and Tolman metrics is intimately linked to the regularity 
conditions of the latter and, for the sake of clarity and completeness of this work, I shall 
discuss briefly the basic argument for Tolman regularity at r = (Bonnor 1974), with 
some additional remarks. Considering a displacement in the 2-surface t = constant, <fi = 
constant, equation ([!]) becomes 

-d£P = ^ (dr 2 + ^fde^ . 

Near r = this 2-surface must be Euclidean, therefore, R' 2 f~ 2 — ► constant 0) and 



f 2 R 2 



R 

that is, 



— ~ r 2 , 



R' 2 r 2 



r 



which is Bonnor's (1974) equation (2.8). Now, if we suppose that / is approximately 
constant for small r, equation (26) can be integrated to get R = r*. Therefore, 

R' 

lim — = limr-^ 1 . (27) 

r^O f r->0 



3 This is especially true as in this problem we are seeking a fractal density profile, and, so, what is 
relevant to the problem are the intermediary results of the integration and not its final values. 
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If / < 1 the limit above becomes infinity. If / > 1 the limit is zero. Hence, R'f 1 can 
only converge to a finite and non-zero constant if / = 1 as r — > 0, and R ~ r also near 



r = 0. So the function f(r) as an exponent in equation (27) and with the value / = 1 
appears to assume a role of criticality in regular Tolman solutions. 

To calculate the matching between flat and Tolman spacetimes, let us start with two 
Tolman regions 

(dS^ 2 = dt 2 - (R\) 2 (f 1 y 2 dr 2 - (Ri) 2 dQ 2 ; r < e; (28) 
(dS 2 ) 2 = dt 2 -(R' 2 ) 2 (f 2 y 2 dr 2 -(R 2 ) 2 dn 2 ; r > e. (29) 

The Darmois junction conditions at the comoving surface r = e were obtained by Bonnor 
& Chamorro (1990) as 

Rr(e,t) = R 2 (e,t); A(e) = / 2 (e); F 1 (e)=F 2 (e). (30) 

If we take region 1 as flat with static boundary, we can write R\ = r, fx = 1, and assuming 



Fi = F 2 = F(r), provided F(0) = 0, the junction conditions (30) are equivalent to 



R 2 [e,t(e)}=e, f 2 {e) = l. (31) 

This result means that on the curved side of the joining surface the solution must be 
parabolic. Considering the parabolic solution ^ the matching conditions (31 ) are equiv- 
alent to solving the following equation for e: 

9F(e) [(3(e) - e} 2 - 8e 3 = 0, F(e) > 0, (32) 

where the geodesic for r < e is the flat one t = —r. Therefore, we can start the numerical 
integration with the initial conditions 

*(e) = -e, 1(e) = 0, (33) 

while 

R(e) = e, R'(e) = 1, R(e) = 0. (34) 

It was the explicit values of the function R[r, t(r)] on the surface r = e that were missing 
in the previous discussion made in paper I. Notice that the junction conditions do not 
require R' to be continuous across the joining surface. 
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It is still possible to simplify the problem by noticing that e = is actually a solution 
of equation (32). Hence, instead of initial values (33) and (34) we have 



t(0) = 0, 7(0) = 0, 



(35) 



and 



77(0) = 0, 77'(0) = 1, 77(0) = 0. 



(36) 

In this way the previous difficulty of indetermination of the function 77 at r = is avoided 
as this method allows us to obtain explicit values of 77 at the origin. This is obviously a 
consequence of the complete formulation and solution of the interface between the two 
spacetimes. 



The initial conditions (35) and (36) may be obtained alternatively without mention 
of the flat-curved junction conditions. Around r = the regularity conditions /(0) = 1, 
F(0) = reduce the field equation ^ to simply 77 = 0. As we know that 77 ~ r when 
r — > 0, we obtain at once the two other values in equations (36). 



In conclusion, we have in practice three sets of initial conditions, where each one, in 



theory, enables the numerical integration of equations (11) and (19): 



(i) the initial values (33) and (34) with e being found by solving numerically equation 



(32); 



(ii) the initial conditions (35) and (36) 



(Hi) assuming initial values (33) and taking e small (~ 10 3 ), but calculating 77(e), 77' (e), 



and 77(e) by putting e directly into equations fel), (pj) and (pi) and their derivatives. 



This third procedure implies neglecting the results (34) of the matching and that will 



almost surely give rise to discontinuities in 77(r, t). Nevertheless, provided e is sufficiently 



small (but not too small otherwise 77 blows up) and equations (11) and (19) are numer- 



ically stable, we would usually get in the end the equivalent results as if we had used 



initial values (35) and (36) with a small perturbation e. As we shall see next, this last 



remark will prove to be most valuable in acquiring confidence in the numerical results. 
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3.2 Numerical Algorithm 

Once the initial value problem for the numerical integration of the first-order ordinary 



differential equations (11) and ( 19 ) is determined, the next step is to choose the numerical 
integrator that propagates the solution over the interval. In this respect it is essential to 
bear in mind that what is being sought is a qualitative fractal behaviour for the solutions, 



that is, solutions which graphically follow the de Vaucouleurs density power law (25). 
Therefore, high accuracy is not important in this problem. In addition, although the 
problem is mathematically tricky, the numerical tasks are relatively simple and, so, 
computational efficiency is of no great concern. In view of this the explicit fourth-order 
Runge-Kutta method with adaptive stepsize control was the integrator chosen (Press et 
al. 1986). 

As discussed in paper I, the aim of modelling a fractal distribution by Tolman's 
solution is to make use of the freedom of the arbitrary functions in order to ascertain 
particular functions f(r), F(r), (3(r) such that the volume density takes the de Vau- 



couleurs density power law (25). The way to check whether a fractal distribution was 



modelled is by taking the logarithm of equation (25) 



log p v = a x + a 2 log d t (37) 

and finding the constants a\ and a 2 by linear fitting over the points obtained through 
numerical integration. The fitting is considered a success if it is linear with a negative 
slope. That can be concluded by visual inspection of the graph and by an acceptable 
goodness of fit. Inasmuch as what we are seeking is a qualitative fractal behaviour, very 
often a visual inspection will be enough. 

I shall describe below in general terms the algorithm used to apply the model and 
obtain the numerical results desired. 

1) START by choosing special forms for the functions f(r), F(r), /3(r), such that 
f(0) = 1 and F(0) = 0; 

2) discretize the radial coordinate rj (i = 1, 2, . . . , n — 1, n; e < r < S ; r\ = e, r n = 
E ; total of n — 1 stepsQ 



In practice the discretization fj will be performed automatically by the adaptive stepsize control 
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3) make i — 1 and choose a set of initial conditions to get t\ and Jj.; 

4) evaluate f h F h (5 h f\, F' h (3',; 

5) if f 2 > 1 or f 2 < 1 then evaluate 6j and consider ti to calculate Ri, R[, R'i, else 
evaluate directly Ri, R[, R'i considering tf, 

6) evaluate z u d h , N Ci , V h p Vi ; 

7) by using the numerical integrator chosen, advance solution t i+ i then use value ti to 
advance solution 

8) make j — i + 1 then i = j; if i = n then continue below; else return to stage 4; 

9) make final evaluations f n , F n , R n , R' n , z n , etc; 

10) with the collection of points p Vk and d\ k (k = 2,3, ... ,n) perform a n — 1 points 
linear fitting according to equation (37); 

11) if fitting is unsuccessful then discard used functions /(r), F(r), f3{r) and go to stage 
1; else STOP. 

The computer routine which performs this algorithm was written in standard Fortran 
77 in double precision and the interested reader can find it in Ribeiro (1992c). Each run 
of the routine, either successful or not, constituted a different numerical simulation. The 
program uses many subroutines of Press et al (1986) with some minor modifications. 
Round-off and truncation errors are monitored at each step by calculating the deviation 
from zero of the field equation (|2j) and its first partial derivative in r. Stability and 
accumulation of errors are monitored by taking the final results as initial conditions to 
reverse the whole integration and compare the results with the initial values. Error prop- 
agation formulas are also used to see whether any of the derived observational quantities 
acquires unacceptably high errors. The units too are chosen bearing in mind the need to 
minimize errors and avoid large numbers. Hence, distance is given in Gpc (10 9 pc) and 



routine which keeps the accuracy of the numerical results within a desired predetermined value (Press 
et al. 1986). 
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in order to have c = G = 1, time is expressed in 3.26 Gyr units (1 yr = 3.16 x 10 7 s) 
and mass is expressed in units of 2.09 x 1O 22 M . 

The first numerical simulations used the set (in) of initial conditions described above 
with e = 1CT 3 . The fractional accuracy for the local truncation errors in the variable- 
step routine was set at 10~ 4 , and this was the value in all simulations, that is, the 
results throughout this paper were all obtained with this accuracy. This set (in) of 
initial conditions applied to equations (11) and (19) led to numerically stable results, 
with the routine taking a relatively small number of steps to go through the interval of 
integration. Around 10 to 20 successful steps were usually enough to cover that interval, 
with 1 or 2 failed steps (see Press et al. 1986, p. 555, for details of the workings of 
an adaptive stepsize control routine). The reverse integration, that is, taking the final 
results to integrate backwards in order to compare with the initial conditions, produced a 
few more failed steps, but the comparison with the initial conditions showed results well 
within the desired accuracy. These details are important to highlight the smoothness 
and stability of the solutions, although this approach actually gave rise to discontinuities 
in R across the surface r = e, as expected (more details about the way the routine works 
and its output are available in Ribeiro 1992c). 

As far as numerical stability is concerned, the situation changed dramatically when 
the set (ii) of initial conditions was used, that is, taking e = and assuming initial values 
(35) and (36). The solutions became unstable about the origin, with the program being 
forced to decrease the stepsize so much that it needed to take on the order of hundreds 
of steps to cover a tiny interval very close to the origin. Once they emerged from this 
unstable region, the solutions started behaving like the previous approach, covering the 
remaining region very quickly. The reverse integration behaved likewise, with hundreds 
of failed steps. Despite this instability the solutions were still very smooth and the 
comparison of the initial conditions with the results of the backward integration were 
still well within the desired accuracy. Thus it would appear that with this approach the 
numerical problem is stiff across the joining surface. 

Integrating a stiff problem with an explicit Runge-Kutta adaptive stepsize control 
routine greatly reduces the efficiency of the code as the initial stepsize chosen puts the 
method at or near numerical instability, which causes a large truncation error estimate. 
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That leads the routine to reduce substantially the stepsize to keep the local truncation 
error within the predetermined accuracy. Then the method is able to successfully inte- 
grate the problem, but uses a far smaller stepsize and a far greater number of steps than 
seems justified if we consider the smoothness of the solution (Lambert 1991). 

Despite the loss of efficiency, our simulations were usually sufficiently straightforward 
and computation time did not pose a barrier in most of them. In addition, double 
precision usually provided enough manoeuvring space for the decrease of stepsize. The 
situation, however, was far worse when the simulations were carried out using the set 

(i) of initial conditions, that is, with e ^ 0. The stiffness was so severe that after a few 
iterations the step length became so small that the difference r i+ i — r, went beyond the 
machine's precision and the routine achieved an apparently zero stepsize. Obviously such 
a case cannot be integrated unless an implicit method is used or, perhaps, if a quadruple 
precision routine is implemented in our explicit Runge-Kutta method. In any case, as 
the set (ii) of initial conditions provides a simpler approach and allows the stiffness to 
be handled by the program in most cases of interest, these were the conditions used to 
obtain all results shown in this paper, unless stated otherwise. 

It is important to mention that when comparing the results obtained by using the set 

(ii) with the ones produced by the set (Hi), the final results of both integrations usually 
differ by about which substantiates the interpretation at the end of §3.1 that the 
results obtained with initial conditions (Hi) are usually equivalent to the ones obtained 
with conditions (ii), but perturbed by about e (see Ribeiro 1992c for an example of such 
a situation). 

Stiffness is associated with the existence of transient terms in the solution: they 
contribute very little to the solution, but the usual methods require that they be ap- 
proximated accurately in order to maintain stability (Gear 1971; Ortega & Poole 1981). 
In this case it appears that the imposition of the regularity conditions upon the Tolman 
solution, or the similar process of joining it to a flat region, give rise to these transient 
terms. Nonetheless, the unanswered question is whether the Tolman stiffness is a mere 
numerical phenomenon in the solutions, with no major implications, or if it indicates a 
deeper physical reason in the Tolman models. It is known that some chaotic dynamical 
systems present stiff behaviour, but the interaction between stiffness and chaos is still 
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unclear (Cartwright & Piro 1992). More work is necessary in order to clarify this point. 

4 The Spatially Homogeneous Special Cases 

In this section I shall present the results of the integration of the spatially homogeneous 
Friedmann models, obtained as special cases of the Tolman metric. The aim of this 
section is twofold: first to extend the analysis of the Einstein-de Sitter model presented 
in paper II to the K — ±1 Friedmann cases, and second to provide numerical results 
of the Einstein-de Sitter case that can be compared with the exact solution already 
presented in paper II. 

As discussed in §2 the Einstein-de Sitter model is obtained as a special case of the 
Tolman solution when / = 1, F = |r 3 , (3 = (3 , where f3 is a constant, and this 
assumption means that along the past light cone the volume density and luminosity 
distance are given by (see paper II) 

_ (3/V /3 - r) 6 
Pv " 54vr(3/? ) 4 ' 

, 9r/3 4 / 3 

0,1 = iy. 

(3/V /3 -r) 2 

Throughout this section I shall assume the value 75 km s _1 Mpc^ 1 for the Hubble 
constant, and this means (3q = 2.7 for the Einstein-de Sitter model in our adopted 
units. Figure 1 shows the logp„ versus log di graph obtained when using the analytical 
expressions above in the parameter range 0.001 < r < 1.5, and figure 2 shows the 
same results obtained by integrating numerically the Einstein-de Sitter model in the 
same parameter range. By comparing the two graphs one can see that the numerical 
integration does reproduce the analytical results. Figures 3 and 4 show respectively the 
null geodesic t(r) and the function I{r) obtained numerically for this case, and there 
the smoothness and almost linearity of the solutions within the range of integration are 
clearly visible. The analytical expressions for these two functions, obtained in paper II, 
may be written as 

*= (/V /3 -0 3 -A), 
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and by means of a few direct numerical substitutions one can see that the numerical 
results shown in figures 3 and 4 do reproduce the analytical functions above, as expected. 

The Einstein-de Sitter model also exhibits stiffness, but this is not visible in figure 
2 because only steps bigger than 0.0025 were stored in the integration shown in that 
figure. Figure 5, however, shows the same integration of figure 2, but with all points 
used by the code being printed. The stiffness is then clearly visible as the program takes 
an excessive number of steps around r = despite the smoothness and linearity of the 
solutions t(r) and I(r) in the range of the Einstein-de Sitter model shown in figures 3 
and 4. 

The open Friedmann model with f2 — 0.2 is obtained by taking / = coshr, F = 
sinh 3 r and f3 = 3.6. Figure 6 shows the numerical results of \ogp v against \ogdi for 
this model, also integrated in the range < r < 1.5. By comparing figure 6 with figure 
2 we can see that in the open Friedmann model p v departs from a constant value at 
about logdi = —0.6 (di « 250 Mpc), while in the Einstein-de Sitter case that happens at 
about log di = —1 (di 100 Mpc). In other words, the departure from a homogeneous 
distribution, where p v is constant, is deeper in the open Friedmann model. Stiffness is 
also present in this case, but in a more severe form as the code took a far greater number 
of steps for the integration to emerge from around the origin than in the equivalent 
Einstein-de Sitter integrations. 

Finally, figure 7 shows the graph for the recollapsing closed Friedmann model with 
fi ^ 4, which is obtained as a special case from the Tolman metric when / = cosr, 
F = sin 3 r, f3 = 0.7. The numerical integrations were also carried out in the same 
parameter range as the other two cases, and figure 7 shows the results of p v vs. d\. In 
this case the volume density starts to change from a constant value at around \og di = 
—1.5 (d t ~ 30 Mpc). Therefore, in the recollapsing model the integration along the 
backward null cone reaches different and earlier spatial sections at closer distances than 
in the other two spatially homogeneous models. Notice that the values of the luminosity 
distance where the departure from homogeneity begins are dependent on the Hubble 
constant used in these integrations, and, therefore, the numbers given above should be 
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considered as rough approximations. 

The recollapsing Friedmann model could not be integrated with the set (ii) of initial 
conditions because the stiffness led to very small first steps, where F w and / ~ 1 
when r — > 0. Hence, the two limits of the interval given by equation ^ quickly became 
extremely big numbers and finding by means of equation ^ broke down due to the 
machine's precision limitations. Besides, the numerical instability often led to either 
shell crossing behaviour at very small r or to unacceptably high error accumulation, 
and in consequence this numerical approach could not be relied upon. As a result the 
recollapsing integrations were carried out using the set (Hi) of initial conditions with 
e = KT 3 . 

The results for the K = ±1 Friedmann cases show that they also do not appear to 
remain homogeneous along the backward null cone. In addition, considering the form 
of the graphs of these two cases we can conclude that similarly to the Einstein-de Sitter 
case studied in paper II, the recollapsing and open Friedmann models do not have single 
fractal features along the past light cone as they do not follow the de Vaucouleurs' density 



power law (25). Therefore, considering the IRAS data, all three spatially homogeneous 
Friedmann models appear to have difficulties in modelling those observations as they 
predict homogeneity where it is not observed, and for deeper regions all models start to 
deviate from it. 

The apparent inhomogeneity of the Friedmann models is a consequence of the fact 
that in the present approach densities are expressed along the past light cone, which 
is where observations are actually made. This means that the volume density p v goes 
through hypersurfaces of constant t where each one has different values for the proper 
density. Therefore, although the Friedmann models are spatially homogeneous, the vol- 
ume density will only be homogeneous, that is, have a constant value, in local regions. 
This will only happen while it is being calculated nearby our present constant time 
hypersurface, and once we depart from its neighbourhood the volume density starts to 
change, with the model becoming increasingly apparently inhomogeneous the further the 
calculations go into the past. The change of the volume density along the geodesic is 
due to it being a cumulative density, averaging at bigger and bigger volumes in a manner 
that adds more and more different local densities of each spatial section of the models. 
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Therefore, under this definition the "homogeneity" of the Friedmann models does not 
survive. It is important to realize that this interpretation of the "homogeneous" Fried- 
mann model is closely related to the fractal approach to cosmology adopted in this series 
of papers, inasmuch as the definition for the volume density used here is a direct conse- 
quence of the assumption that the large-scale structure of the universe can be described 
as being a self-similar fractal system. 

As final remarks, when one approaches the singularity the model is usually said 
to become denser because the local density p tends to infinity. However, this is not 
the case for the volume density p v , which vanishes at the big bang. This behaviour 
definitively happens in the Einstein-de Sitter model (see paper II for the asymptotic 
analysis of this model) and due to the similar form of the graphs of log p v vs. log di 
in the open and closed Friedmann models (see figures 6 and 7) as compared to the 
Einstein-de Sitter case (figures 1 and 2), we can then assume that this vanishing volume 
density is a feature of all Friedmann models when they approach the singularity. Such a 



behaviour for p v is, once again, explained by its adopted definition [eq. (23)]: at the big 
bang singularity hypersurface the observed volume (as defined through the luminosity 
distance) is infinite, but the total mass within it is finite. It is interesting to note that 
Wertz (1970) had already postulated a vanishing global density as a requirement for 
the so-called Pure Hierarchical Models, and similarly Pietronero (1987) conjectured the 
same sort of result for a fractal distribution. Therefore, it does seem valid to say that at 
least some aspects of a hierarchical (fractal) cosmology can already be found within the 
standard Friedmannian cosmology, provided one starts with the appropriate definitions. 

Empirically, the effect of zero global density in these models as they approach the 
big bang can be seen as the consequence of the fact that in this work the galaxies one 
sees empirically on the past light cone are plotted at their luminosity distances, and 
the luminosity distances go to infinity as z goes to infinity while the number of galaxies 
remains finite. Observers do often plot their data in ways like that, and an example 
is Geller & Huchra's (1989) slices where galaxies are plotted at a distance proportional 
to their redshift z. If they had extended their plot to z equals infinity only a finite 
number of galaxies would be inside the horizon and one would see the galaxies thin out 
to zero. It is obvious, however, that if a different definition of distance had been used (say 
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comoving distance), this effect of zero global density at the big bang singularity would 
no longer be necessarily found, nevertheless the redshift or the luminosity distance are 
observable quantities, independent of cosmology, and observers do use them. Therefore, 
we can then consider the de Vaucouleurs' density power law as being essentially an 
empirically observed relation when observables are plotted, and from this perspective 
this work can be viewed as an attempt to see whether there are any inhomogeneous big 
bang cosmologies that would allow this empirical relation to be extended to larger scales. 



5 Tolman Fractal Solutions 

Before discussing the solutions themselves, let us first make some remarks about the 
three arbitrary functions which may help to simplify the simulations. The function F(r) 
cannot be a constant, otherwise due to the regularity condition F(0) =0 this constant 
would have to be zero and there would not be any dust in the model. In addition, due to 



equation (10) we have to restrict ourselves to F' > 0, otherwise we would get negative 
mass in the model. Hence F(r) must be a positive and monotonically increasing function 
for r > 0. Notice, however, that there could be ranges for r where F(r) is a nonzero 
constant and in those ranges F' = 0, implying the existence of vacuum regions in the 
model. The only restriction for the function f3(r) is that its values must be such that 
the model is kept within its physical region, that is, such that t + (3 > 0. The aim of 
the simulations is to find a qualitative fractal behaviour which could be expressed in 
terms of a relatively simple metric. Therefore, it is desirable to obtain fractal behaviour 
in terms of simple functions. So, wherever possible, I shall restrict the choice of f(r) 
to the forms it takes in the spatially homogeneous cases, that is, specializing only to 
/ = cosr, 1, coshr. 

In deciding which solutions are acceptable as fractals, we have the mathematical 



criterion of verifying whether they follow the equation (25). However, some basic ob- 
servational constraints must also be obeyed in order to have physically realistic models. 
The first and most obvious is the observed linearity of the redshift-distance relation for 
z ~ 1. The second is the value of the Hubble constant itself which, considering the 
present uncertainty in its measurements, is accepted to be in the range 40 km/s/Mpc < 
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H < 100 km/s/Mpc (see White 1990 for a recent account on the subject). Finally we 
have a constraint in the fractal dimension of the distribution. Considering Pietronero's 
(1987) analysis based on the measurements of the 2-point correlation function, we would 
have D w 1.4. However, as other authors claim different and higher values for D (Deng, 
Wen & Liu 1988; Calzetti & Giavalisco 1991), I shall loosely assume here that the wider 
range 1 < D < 2 is reasonable enough for the purposes of this work. 



5.1 Parabolic Solutions 

The parabolic class of Tolman models is the natural one to start the search for fractal 
solutions, inasmuch as one of the three functions is already specified, that is, / = 1 in 
this class. Nevertheless, and despite the simplifying remarks above, it was not at all an 
easy task to find those results. It took over 250 different simulations to find suitable 
solutions of parabolic type, and one of the major problems was to find the range of the 
constants in the functions where the overall solution does have a power law behaviour 
within the interval of integration. 

Despite those difficulties, fractal parabolic solutions along the past light cone do 
exist, although I could only find them by assuming a non-simultaneous big bang, that 
is, I could not find fractal results with constant (3{r). Below are the particular forms of 
the remaining arbitrary functions that led to fractal behaviour in parabolic models: 

F = arP, , . 

38 

, P = Po + Vo r\ 
and 

' F = ar p , , 

f\ \ 39 

P = In (eP° + rjir) , 

where a, p, q, /3 , r/o 5 Vi are positive constants, [j The simulations showed that in both 
cases above a must be around 10~ 4 to 10~ 5 and p and (3q can vary from around 0.5 to 4. 



For functions (38) q lies around 0.65 and t]q around 50. For solutions (39) r\\ varies from 



about 1000 to 1300. It is difficult to attribute specific roles for each constant as usually 



5 Actually, by a coordinate transformation F can be made equal to |r 3 , so the only genuine arbitrary 
function is (3. 
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a change in one of them changes the whole behaviour of the solution, but we can very 



roughly say that a, q, rj and r)i are more related to the linearity of equation (37) and 
that p and (3q are more related to its slope. In other words, changing p and /3q affects 
more sharply the fractal dimension of the solutions, but without changing too much their 

In 



linearity. In addition, if in equations (39) we use instead f3 
linearity of the solutions is slightly improved. 



e /3o _|_ 7 y 1 _|_ r 2^ 



the 



Figure 8 presents the results of an integration obtained with functions (38). The 
fractal dimension is D = 1.4 and the integration was carried out in the interval < r < 2. 



In figure 9 one can see the numerical results for the functions (39) integrated in the range 
< r < 5, with the resulting fractal dimension D = 1.7. Both figures also show the 



straight lines fitted according to equation (37) and the qualitative power law nature of 
the solutions in the integrated interval is visually obvious. We can also see from the 
figures that the two solutions actually are not much different from one another. 

These solutions, however, fail to meet the other observational criteria outlined above. 



Figure 10 shows the di vs. z plot of results of the integrations of functions (38) and the 
non-linearity of the relation is clearly visible. In addition, the associated Hubble constant 
is far too low and these reasons seem to rule out as observationally unrealistic a model 



based on the functions (38). Figure 11 presents the same results for functions (39) and 
one can see that the redshift-distance relation is mostly linear, especially for the region 
z > 0.011. The Hubble constant associated with the linear part of the diagram is around 
13, far below the lower limit for Hq given above, and considering this point it seems that 



the functions (39) also do not seem to produce a realistic fractal model. 

Two remarks must be made about the fractal parabolic solutions above. In the first 
place, from a mathematical point of view the linearity of the Hubble law ffl must be 
approximately valid through the origin, obeying the equation di = cz/Hq. This means 
the redshift distance relation in figure 11 would not be acceptable even though a section 
of the curve looks linear. However, from a physical viewpoint one might argue that 
the linearity of the Hubble law is not really needed very near the origin where it would 
anyway be messed up by local peculiar velocities. In the second place, although it might 



6 Here by Hubble law I mean the redshift-distance law. See Harrison (1993) for a discussion about 
the distinction between the Hubble law and the velocity-distance law. 
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be argued that a different and more observationally acceptable value for the Hubble 
constant could be obtained by changing the constants until the desired value for H Q 
is produced, it must be stressed that the fractal behaviour of the solutions is sensitive 
(sometimes very sensitive) to the choice of these same constants. In other words, the 
fractal dimension itself and the limits over which the de Vaucouleurs density power law 
applies are sensitive to the choice of these constants. Therefore, changing them in order 
to get a desired value for the Hubble constant will only be succesful at the very likely 
cost of destroying the fractal feature of the solutions over the desired observational range. 
Due to this reason such a procedure was ruled out in all simulations. Q 



As final remarks, although one cannot say that the functions (J38J) and (39) are the 
only ones to produce fractal parabolic solutions as there may be others which lead to 
fractal behaviour but were not considered here, the experience with all simulations makes 
it hard to see how other solutions of this kind could be produced in terms of simple 
functions. This point added to the unrealistic distance-redshift relations of the models 
leads us to conclude that parabolic type Tolman models do not appear to be able to 
produce physically realistic fractal models in the backward null cone. 

5.2 Elliptic Solutions 

The experience with parabolic models was valuable guidance for further studies of Tol- 
man fractal solutions. Elliptic models with power law behaviour in the density were 
found directly from their parabolic counterparts: 

cos r, 

(40) 



and 




/ = cos r, 

F = ar p , (41) 
/3 = In (e* + r/ir) , 



7 This sensitivity of the fractal solutions to small changes in the constants seems to indicate that 
the fractal solutions presented in this paper appear to be structurally fragile (Coley & Tavakol 1992). 
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where a, p, q, r/o, 771 are positive constants which must lie approximately in the same 
range as in the parabolic cases in order to produce fractal solutions,^] apart from a which 
must be around the value 10. We have to remember that in the elliptic class models the 
mass is higher than in the other two classes, which means that the dust does not escape 
from its own gravitational field and this leads to the eventual halt in the expansion when 
every part of the model starts to contract. Therefore, it is reasonable to expect a higher 
value for a in these cases. 

The numerical integrations of all elliptic models were carried out using the set (Hi) 
of initial conditions as described in §3.1. This was done for the same reason as in the 
recollapsing Friedmann model: the set (ii) of initial conditions and equation (JsJ> broke 
down due to severe numerical instability at very small r, which produced unreliable 
results, if any. 



Figure 12 shows the di vs. p v results for functions (40), with a resulting fractal 



dimension of D = 1.7, and figure 13 shows the same results for functions (41) where 
D = 2 was the dimension found. Figure 14 shows the distance-redshift relation for the 



model with equations (40) and we can see that the plot is not linear. Finally, figure 15 
shows the di vs. p v plot for the model of equations (41), and although the relation is 
more or less linear it leads to the low value of roughly 36 for Hq. Considering all those 



results together it also seems that equations (40) and (41) produce models which are not 
compatible with observations. Models with constant (3 have less linear log p v vs. log di 
plots than the two previous ones and lead to too high values for Hq and D, making them 
also incompatible with observations. 



5.3 Hyperbolic Solutions 

The best fractal results in terms of power law behaviour for p v vs. di and agreement 
with observations were obtained by numerical simulations of hyperbolic type solutions. 
The specializations for F(r) and /3{r) studied so far also lead to power law behaviour in 
the density once the function f(r) is modified accordingly. Thus, we have the following 

8 Since one can make a coordinate transformation in r, only two of the functions /, F, (3 are really 
arbitrary. 
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sets of functions: 



/ = cosh r, 

F = arP, (42) 
^ P = p + r] r«, 
and 

/ = coshr, 

F = arP, (43) 
P = In (e* + r/ir) , 

where, as before, a, p, q, Po, r] Q , rji are positive constants. In this case the constants 
will assume approximately the same numerical ranges as in the parabolic case in order 
to produce fractal solutions, apart from p which for values smaller than 1.4 and greater 
than 2.5 makes the problem too stiff to be integrated with initial conditions (ii). For 
those values the use of the initial conditions (in) will produce the desired results. 

Figure 16 shows the p v vs. d\ numerical results of the integration of the model using 
equations (42) and the fitted straight line, which leads to a fractal dimension D = 1.4, 
while figure 17 shows the redshift-distance diagram. Although the latter relation is 
approximately linear, it leads to an Ho pa 23 km/s/Mpc, a value too low to make 
equations (42) a model compatible with observations. 

The power law behaviour, leading to fractal dimension D = 1.3, of the results of 



the numerical integrations of the model formed by functions (43) can be seen in figure 
18. Figure 19 shows the redshift-distance diagram of the same model and we can see 
that the plot is linearly well approximated. More important, a simple calculation shows 
that the slope of the approximated line formed by the points in figure 19 gives Ho = 61 



km/s/Mpc. Therefore, the model formed by functions (43 ) is the first so far to reasonably 
agree with all basic requirements outlined above in order to make the model compatible 
with current observations: linearity of the d% vs. z diagram, Hubble constant within 
the currently accepted uncertainty, power law behaviour for the p v vs. d\ results and 
fractal dimension around the values which are in agreement with present calculations of 
the 2-point correlation function. Note incidentally, that the integrations were stopped 



As in the elliptic case, here again only two functions are really arbitrary. 
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li = 350 Mpc, which corresponds to z = 0.07, as this is the redshift of the deepest 
self-similar structures identified in the IRAS survey (Saunders et al. 1991). 

Even better results were obtained with the mathematically and physically simpler 
model below: 

/ = coshr, 

F = arP, (44) 
/? = A)- 

This is obviously a simultaneous big bang model, like its Friedmann counterpart, and is 
physically and mathematically the closest model to the open Friedmann one, differing 
only in the function F(r) which gives the distribution of dust. Figure 20 shows the power 



law behaviour of p v vs. di of the model formed by functions (44), with a fractal dimension 
of D = 1.4. Figure 21 is the redshift-distance diagram of the same model where we can 



see the good linear approximation given by functions (44). The slope of the points gives 
Hq = 80 km/s/Mpc, and it is interesting to note that recent measurements made by 
two different methods suggest a Hubble constant very close to this value (Peacock 1991). 
Actually, for f3 = 3.6, the value used to get the results shown in figures 20 and 21, we 
would have an age of the universe of about 12 Gyr, which is a lower limit if we consider 
the age of globular clusters (Peacock 1991). Therefore, also in this point of the age of 
the universe the model (44) agrees reasonably well with observations. The integrations 
with functions (44) were also stopped at z = 0.07, which in this case corresponds to the 
luminosity distance di = 270 Mpc. Finally, figure 22 shows the results for cumulative 
number counting vs. redshift produced by the model under consideration. 



6 Discussion 



6.1 A Metric for a Smoothed-out and Averaged Fractal 

In the previous section it was shown that Tolman fractal solutions do exist, and that the 
only ones compatible with observations are the hyperbolic type solutions obtained by 
means of the specializations given by equations (43) and (44)]^] As in equations (43) the 



10 Notice however that fractal solutions of elliptic and parabolic types may, in principle, be obtained 
from other more complex specializations of the arbitrary functions than the ones considered in this 
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function (3{r) is not constant, this means that the model has no simultaneous big bang. 
In other words, in a model of this sort the big bang singularity hypersurface occurred at 
different proper times in different locations, and the age of the universe is different for 
different observers at different radial coordinates. More specifically, as /3(r) in equations 
(43) is an increasing function, regions at smaller r are younger than at bigger r, and the 
youngest region of the model is "here", at r = 0. 

An universe model where some regions are older than others is not as odd in terms of 
accepted ideas of galaxy formation as it might seem at first. Inasmuch as the observed 
universe is lumpy, in a spatially homogeneous Friedmann universe where p = p(t) and 
the big bang singularity is simultaneous, there must be density fluctuations Sp/p of 
some kind in order to form galaxies, and it is necessary to have some sort of metric 
perturbations for that to happen. So, at the era of galaxy formation, which may be 
defined as a hypersurface of constant time in order to agree with our intuition in an 
unperturbed metric, in the perturbed metric the overdensities p Q occur at time t D and 
the underdensities p u occur at t u , and t Q must be different from t u , otherwise there 
would not be any fluctuation at all. In other words, due to the density fluctuations 
Sp/p, in the perturbed metric a hypersurface of constant density no longer coincides 
with a hypersurface of constant time. Therefore, a deviation of the Friedmann metric 
from spatial homogeneity, even if it is small, is essential for lumpiness and, hence, some 
regions will inevitably have different local times than others. In the perturbed metric we 
may even define the era of galaxy formation as being a hypersurface of constant density. 
In other words, a non-simultaneous big bang seems inevitable in order to form galaxies 
in the standard scenario, even if those differences in local times are small. 

Note that this discussion assumes that a hypersurface of simultaneity is defined by a 
specific value of the proper time, which is a logical thing to do in an unperturbed metric. 
In a perturbed metric, however, one could define the hypersurfaces of simultaneity in a 
different way, which would mean a different choice of gauge by which the perturbed and 
the non-perturbed spacetimes are related. 

Nevertheless, considering it is desirable that fractal models be as close as possible 
to their Friedmann counterparts, and also considering mathematical simplicity, I shall 



paper, and these solutions could be compatible with observations. 
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take the specializations given by equations (44) as the best modelling of a relativistic 
hierarchical (fractal) cosmology by Tolman's spacetime. Let us now write this model 
explicitly. Its metric is expressed as 

dS 2 = dt 2 - ( R 2 ) dr 2 - R 2 (d6 2 + sin 2 ## 2 ); r > 0, R[r,t(r)] > 0. (45) 
\cosh r J 

The Einstein field equation for this metric may be written as an energy equation (Bondi 
1947; paper I) 

*-U{R)=E{r), (46) 



where 



is the effective potential energy, and 



U(R) = ^ (47) 



E(r) = -sinh 2 r (48) 



is the total energy within r. The solution of equation (46) is 



ar 



p 



(cosh 26-1), (49) 
8F(r) 

where 6 is given by 

4 [t(r) + fa] [2,E(r)] 3/2 = ar p (sinh 26 - 26) , (50) 



t(r) is the solution of the past radial null geodesic 

dt R' 
dr cosh r ' 

and the local density is expressed as 

n2 



(51) 



p = (52) 

7rarP+ 1 J R / (cosh26 - l) 2 

The essentially new physical feature of the model above is its single difference from the 

open Friedmann one: the form of the function for the gravitational mass. That is given 

by F = ar p , while in Friedmann this function must be F = 6 2 sinh 3 r. Remembering 

that a = 10~ 4 — 1CT 5 and p — 1 — 2.5, the fractal metric (45) appears to have a more 

rarefied dust than its Friedmann equivalent. As fractal models are characterized by a 

power law nature of their average densities, with fractional exponents smaller than 3, it 

is hardly surprising that such models have a more rarefied distribution of mass. 
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6.2 Evolution of the Fractal Model and Comparison with the 
Spatially Homogeneous Case 



The constant (3q in equation (50) gives the universal big bang time if we define "now" 



as t — 0, and this fact allows us to study the evolution of the fractal model (45) through 
the easily computable manner of simply varying /3 - This investigation permits us to 



answer the question of whether or not the fractal features of the metric (45) are present 
at different epochs. Bearing this point in mind, I carried out simulations for different 
values of (3q in the interval < r < 0.07, but keeping a = 10~ 4 and p = 1.4 in all of 
them. I found out that the model under consideration remains fractal in integrations 
where (3q = 1.5, 2, 2.5, 3.6, 4.5 and 6, with the resulting fractal dimensions being D = 1.5, 



1.5, 1.5, 1.4, 1.4 and 1.4, respectively. These results show that the metric (45) effectively 
models a fractal distribution of dust at different epochs, with a remarkable constancy 
in D. All those integrations had r = 0.07, the final value of the integrating parameter, 
corresponding to z = 0.07, although the luminosity distance varied from d\ = 110 Mpc 
when /3q = 1.5 to di = 450 Mpc when f3o = 6. This variation may be physically explained 
as due to changes in the Hubble constant itself, whose values get bigger at earlier epochs. 
Finally, those results together with some simulations on different values of p suggest it 
is possible to propose a simple, but very restricted, relationship between p and D. For 
< r < 0.07, 1.5 < (5q < 6 and 1.4 < p < 2.5 we can say approximately that D = p±0.1. 



An interesting question about the model (45) is to see how it would compare with 
the spatially homogeneous open Friedmann one. Looking at the results of the latter in 
figure 6 and the former in figure 20 we can qualitatively see that the absolute value of 
the difference in log p v between the two models starts as zero, but increases rapidly. In 
the analytical study of the Einstein-de Sitter model presented in paper II it was shown 
that at the big bang singularity hypersurface the volume density p v vanishes and the 
luminosity distance di goes to infinity, and this effect is a consequence of the definition 



(23) of the volume density: at the big bang the volume (22) is infinity, but the total mass 
is finite. As already said at the end of §4, we can therefore expect a similar asymptotic 



effect in both the open Friedmann case and the fractal model (45), and this means that 
the difference in p v between these two models should start decreasing after reaching a 
maximum in its increase. Actually, figure 6 already shows a sharp decrease of p v in 
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the open Friedmann case, but as the integrations of the model (45) in figure 20 did 
not go far enough in the past, the results cannot show where this maximum might be. 
Nevertheless, based on this reasoning we can deduce that the observational relations of 



the fractal model (45 ) appear to be asymptotically Friedmann when calculated along the 
past light cone. 



6.3 Relations to Homothetic Self- Similarity 

Another topic which deserves some investigation is the relation, if any, between self- 
similarity due to fractals and self-similarity due to homothetic Killing vectors (Cahill & 
Taub 1971). There is some interest in this point because some attempts have been made 
to explain large-scale voids and clusters by self-similar perturbations of a Friedmann 
universe (Carr & Yahil 1990), and also because the first attempt to propose a workable 
relativistic hierarchical cosmology was done assuming a homothetic self-similar metric 
(Wesson 1978, 1979). 

In general relativity a spacetime is called self-similar if all metric components can be 
put in a form in which they are functions of a single independent dimensionless variable 
which is a combination of the spacetime coordinates. Mathematically, this corresponds 
to the existence of homothetic Killing vectors, meaning that spherically symmetric simi- 
larity solutions are unchanged by coordinate transformations of the form t — > bt, r — > br, 
for any constant b obeying the conformal transformation g^r, t) —>■ p-g^r, t) (Cahill & 
Taub 1971). Physically, spherically symmetric similarity spacetimes with A = contain 
no fundamental scales or dimensional constraints (Henriksen & Wesson 1978), and it 
would seem that these features make homothetic self-similarity a possible mathematical 
version of the scaling idea behind the empirical hierarchical clustering concept. However, 
homotheties, as defined in general relativitity, are basically geometrical features which 
will not necessarily translate themselves in observable quantities. In other words, the 
geometrical scaling features of the model do not necessarily mean that its observational 
relations are also scaling. For this reason it seems that homothetic self-similarity provides 
an unsatisfactory manner of modelling hierarchy. 

It is beyond the aims of this work to make a general discussion about the possible 
relations between these two types of self-similarity. However it is of interest to investigate 



29 



whether or not the fractal solutions presented in §5 are homothetic. Recently Lemos & 
Lynden-Bell (1989) and Ponce de Leon (1991) studied homotheties in Tolman models and 
found that specific criteria are necessary to maintain the assumed similarity symmetry 
in the solutions. For models where A = (our case here), they showed that the first 
criterion is for the function f(r) to be constant, that is, each comoving shell must have 
the same total energy. That immediately tells us that all fractal solutions of elliptic 



and hyperbolic type studied here, including the metric (45), are not homothetic, leaving 
only the parabolic solutions still to investigate. The next criterion says that homothetic 
solutions with f(r) = 1 restrict the mass distribution to have the form F(r) = constant x 
r 2rj+3 , where p is a constant, and this is the case in the solutions (38) and (39). The 



final requirement for homothetic symmetry to hold in our Tolman solutions demands 
that for / = 1 and p ^ the big bang hypersurface must be of the form /3(r) = a + br~ p , 
where a and b are constants. The solution ([381 satisfies these three requirements only 



if p = 3 — 2q. The solution (39) can have (3{r) reduced to its form in equations (38), 



but then the value p = 1 must hold to satisfy the third requirement. Therefore, very 



restricted cases of the fractal solutions (38) and (39) have also homothetic self-similarity. 

In conclusion, from what we have seen it does appear valid to say that fractal self- 
similarity is a much weaker requirement on the solutions than homotheties, and although 
we have reached this conclusion looking in more detail only at the Tolman spacetime, 
based on the self-similar requirements for both cases it seems reasonable to suppose that 
this conclusion may well be valid in general. 

6.4 Fitting Condition 

In the previous section of this paper we dealt with the problem of finding a specific Tol- 
man model which best represents the observed inhomogeneous distribution of galaxies, 



and in that respect it was concluded that the metric (45) is the simplest one to achieve 
this aim. In other words, what was being sought was the optimal way of fitting the 
Tolman metric to the real lumpy large-scale structure of the universe. It was discussed 
in paper I that it may be desirable for us to have a Friedmann metric as background 
spacetime to the inhomogeneous Tolman region (or regions) used here to describe a frac- 
tal distribution of galaxies, and having found the specific forms for this inhomogeneous 
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region an important question arises at once: what are the implications that the fractal 



metric (45) and the hyperbolic solution (43) bring to a possible Friedmann spacetime 
background? Answering this question is equivalent to finding a response to the "fitting 
problem" in cosmology, in the specific context of this work. 

Ellis & Stoeger (1987) have outlined the fitting problem as being the search for an 
ideal Friedmann model which best fits another cosmological model that gives a realistic 
representation of the universe, including all inhomogeneities down to some specified 
length scale. In Ellis & Stoeger's words, "the approach resembles that used in geodesy, 
where a perfect sphere is fitted to the pear-shaped earth; deviations of the real earth 
from the idealised model can then be measured and characterized". Various ways in 
which this approach might be tried are discussed in detail by them, however here we 
shall restrict ourselves to the specific one outlined by Ellis & Jaklitsch (1989) where 
the matching between the Tolman and Friedmann spacetimes is interpreted as a fitting 
condition. 

It was shown in paper I that the Darmois junction conditions between the two space- 
times under consideration require that / = g' on the joining surface. Here g = sinr, r, 
sinh r is the function which determines the curvature of the Friedmann model. As both 



the solution (43) and the metric (45) are of hyperbolic type with / = coshr, there is 
no way of satisfying this condition when r > unless we have g = sinhr. Therefore, 
the first response to the fitting problem in this context says that our Tolman fractal 
solutions imply an open Friedmann background model. 

It was also shown in paper I that the matching between these two spacetimes severely 
restricts the gravitational mass inside the Tolman cavity. If m(r) = F(r)/4 is the 
gravitational mass of the Tolman region within a comoving radius r, and m(x) = 
47r/ia 3 (t)g 3 (x)/3 is its Friedmann equivalent for a radius x and dust density /x, the junc- 
tion conditions demand 

m(S ) = m(S ), (53) 
where S is the constant that defines the joining surface r = x = S between the two 



spacetimes. Therefore, as discussed by Ellis & Jaklitsch (1989), equation (53) allows us 
to choose the Friedmann background model whose density is appropriate to our lumpy 
Tolman model. Let us see in more detail how this background spacetime can be specified. 
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In the open Friedmann model the local density is given by 

3b 2 



167ra 3 (f)' 



and as F = ar p in our fractal models, we can then write the equation (53) as 



aE p = 6 2 sinh d S . (54) 

The value of the parameter b 2 is what we are seeking in order to determine precisely 
the Friedmann background and give a more accurate answer to the fitting problem in 



this context. Equation (54) shows that b 2 is dependent on the other three parameters 
a, p, Eq, and hence, there is a certain degree of flexibility in choosing the mass of the 
background spacetime. Thus, even when the interior region is determined by known 
values of a and p, different values of 62 are obtained according to exactly where the 
joining surface is located. 

We can work out how the Friedmann background is in the case of the numerical 



integrations of the model (44) shown in figure 20. We have a = 10~ 4 , p = 1.4, and if 
we take E = 0.07 (the value where the numerical evaluation ends) we get b 2 = 0.007, 
Q = 0.002 and Hq = 83 km/s/Mpc. This low value for Qq is in the lower limits of the 
interval where it has been reportedly measured. However, it is important to notice that 
in the approach used in this work no kind of dark matter was considered, but only the 
luminous matter associated with galaxies. Galactic luminous matter gives a value for 
of the same order of magnitude as the one found for the background model above (see 
White 1990, p. 38). 

As a final remark, we have so far considered the interior Tolman region joining di- 
rectly to the exterior Friedmann metric. That does not need to be always the case and 
we can envisage an interior region surrounded by one or more intermediary regions with 
higher or lower densities, in a scheme designed to model specific observational features. 
For example, a structure like the "Great Wall" (Geller & Huchra 1989; Ramella, Geller 
Sz Huchra 1992) could be modelled by an intermediary overdensity region before the 
background spacetime is reached, and with an underdensity interior fractal Tolman re- 
gion (see Bonnor & Chamorro 1991 on how to join an underdensity Tolman region to an 
overdensity intermediary section). Such a modelling will obviously increase the value of 
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fio for the fitted Friedmann background. This method, however, demands more detailed 
work in order to achieve a model where this sort of structure is precisely characterized. 
I shall not pursue further this study here. 

6.5 Criticism of Criticisms of Fractal Cosmology 

As the final issue of this section, I shall discuss some of the objections raised by some 
authors against a fractal cosmology. The first type of criticism is contrary to a possible 
unlimited fractal pattern for the large-scale clustering of galaxies, and although some of 
the critical voices do accept fractals at small scales, their objections are usually based 
either on reasoning from the 2-point angular correlation function (Peebles 1989), or 
on supposed strong theoretical limitations of the standard Friedmannian cosmology. 
Therefore, in one way or another those authors see the strong need for a crossover to 
homogeneity on the fractal structure, at a scale yet to be agreed upon. 

Criticisms based on the angular correlation function have been addressed by Coleman 
& Pietronero (1992) and I shall not discuss them here, although this point was briefly 
mentioned in paper II. The theoretical criticisms, on the other hand, must be addressed 
in this paper, and for this purpose I shall reproduce here two quotations from Martinez 
(1991) which well represent this point of view, although he is by no means the only one 
to raise such kind of objections. Martinez states that "... it should be noted that in 
the standard cosmology, the distribution of mass must tend to a non-zero finite density 
when averaged over large volumes" [^j From a relativistic point of view, the problem 
with this statement is its failure to specify where this average is supposed to be carried 
out. It is correct to say that if in a Friedmannian cosmology we make averages of 
density at spacelike hypersurfaces of constant t, those averages cannot be zero as this 
cosmology is spatially homogeneous. However, at large volumes such averages would 
be observationally irrelevant as astronomy in the electromagnetic spectrum is actually 
made along the backward null cone and not at such spacelike hypersurfaces. Therefore, 

11 The expression "large volumes" used in this quotation is imprecise, and can be interpreted as 
meaning cither big local volumes or volumes which are big enough to be no longer considered as local. 
In a cosmological context the latter is more appropriate and from now on I shall assume the expression 
"large volumes" to mean non-local ones. 
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the statement above is only true in an unobservable situation, and considering that voids 
and clusters of galaxies were, and still are being, identified in the so-called redshift space, 
which lies on the past null cone, this is where the average of density must be carried 
out. Hence, the quotation above is inappropriate as an objection to a fractal structure 
for the distribution of galaxies. These two types of averages will only coincide locally, 
and it will depend on the model and the value assumed for the Hubble constant in 
order to establish what scales are local, although, in any case they will certainly differ 
at large volumes. In effect, in this imprecise formulation this statement may actually 
reinforce, or be taken by, a common misinterpretation of the standard model, which, due 
to inappropriate Newtonian analogies, confuses the model's geometry with its observable 
quantities . 

Martinez (1991) goes on and states that "... a fractal universe without a crossover 
to homogeneity (...) implies a vanishing density for very large volumes and this idea 
cannot be accepted without creating important additional problems". It was shown in 
paper II (and rediscussed in §4) that if averages on density are made along the backward 
null cone, at the big bang singularity hypersurface the luminosity distance goes infinite 
and this average density is zero. That happens in the Einstein-de Sitter model, the most 
popular of the cosmological models, and no additional hypothesis or change in the metric 
was done to achieve this result. The point is, once more, where the average is made and 
which definition of density is adopted. Thus, again the statement is in fact an untenable 
objection since the standard cosmology does have a vanishing average density without 
any important additional problem. Having or not having zero average density in the 
model is just a question of interpretation. 

It should be clearly understood that the above criticisms of Martinez's (1991) state- 
ments are made solely on the grounds of the standard Friedmannian cosmology, and 
there is no need whatsoever to mention any fractal hypothesis in order to show the im- 
preciseness and inappropriateness of such statements. The important point being that 
even accepting the spatially homogeneous standard Friedmannian cosmology, this model 
tells us we would only be able to see, through our telescopes, its homogeneity locally. 

Closely related to this point of local homogeneity is the issue, one could argue, of 
how we would understand in this context the reported uniform distribution of some 
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deep samples like radio sources. In the first place it must be made very clear that the 
discussion made so far is aimed at showing that from a theoretical point of view there 
is no constraint to an unlimited fractal distribution, even from within a Friedmannian 
framework, but that does not imply the fractal system is indeed limitless, as an upper 



cutoff to homogeneity is not yet ruled out. Nonetheless, a simple calculation in the 



Einstein-de Sitter model, as presented in paper II, will show that a 30% decrease in p v 
(from the value at present time t = 0) occurs at z ~ 0.1 (d; ~ 500 Mpc), and this 
means that even considering such high error in the determination of the volume density, 
this is, roughly speaking, the maximum approximate range where the homogeneity of 
this model could be observed. Beyond this range the homogeneity of the Einstein-de 
Sitter model would no longer be observed in the past light cone. Thus the first issue 
raised by this result is a problem of methodology: curvature effects occur in Friedmann 
models at much closer ranges than usually assumed, and this means that those surveys 
must consider in their data analysis expressions along the past light cone. Currently, 
calculations of observational relations where the backward null cone is taken into account 
is a very much neglected problem in cosmology. Secondly, if it is confirmed that in those 
deep surveys the distribution is really uniform, that would put the Friedmann model in 
even greater difficulties as it would appear to predict inhomogeneity in deeper ranges 
where this would not be observed. Thirdly, the sceptical viewpoint on this issue would 
be to argue that usually deeper observations are less precise than shallower ones, and 
previous claims of the so-called homogeneous "fair-sample" finally being observed did 
not stand once more refined and complete observations were made. Historically, the 
range at where the homogeneity is, or would be, finally reached has being pushed further 
and further away as more complete data become available and observational techniques 
improve, and so, the sceptic may say, we may not have necessarily seen the end of this 
story. 

In addition to the points discussed above, a second kind of criticism to the fractal 
cosmology has been voiced by Peebles, Schramm, Turner & Kron (1991) in the following 

12 The proposal for an upper cutoff to homogeneity in the fractal system appears to have been 
initially advanced by Ruffini, Song & Taraglio (1988), although Pietronero (1987) had already made a 
discussion about this issue. 
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form. "If the galaxy distribution had been observed to follow a pure scale-invariant frac- 
tal, (...) the closely thermal spectrum and isotropy of the cosmic background radiation 
in this highly inhomogeneous Universe would have been a deep puzzle". First of all, 
it should be said that in an inhomogeneous model with a Friedmann background, the 
apparent discrepancy between the inhomogeneity of the model and the isotropy of the 
microwave background is not really an issue as the junction conditions already require 
that an overdensity must be compensated by an underdensity before the uniform region 
is reached, in order that the average densities will be the same (see paper I, §4). Never- 
theless, the most important point is the result already obtained in paper II and extended 
in this paper for the other Friedmann models: the standard Friedmannian cosmology 
may be taken to be inhomogeneous depending on how we look at it. That means that 
the apparent contradiction between inhomogeneity and the isotropy of the microwave 
background may not be a contradiction at all. This is an essential point in order to put 
Peebles, Schramm, Turner & Kron's statement above into perspective, as we have already 
seen in this paper and in paper II that at relatively modest luminosity distances and 
redshifts, even the standard spatially homogeneous Friedmannian cosmology becomes 
highly inhomogeneous on the past light cone because p v departs considerably from its 
constant value in our constant time hypersurface the further into the past we look, and 
p is dependent on r. Considering that so far even the deepest all-sky redshift surveys 
have failed to reach the so-called "fair sample" where the homogeneity is supposed to 
be, waiting for us to discover it, perhaps it is about time to ask if the cosmic background 
radiation could be accommodated in a different cosmology, or wonder whether it is really 
a deep puzzle. 

After all this discussion, we are left with an important point to consider. If the 
supposed theoretical need for a crossover to homogeneity in the fractal system is much 
weaker than previously thought, we have the question: is it really necessary? Since 
we have seen here and in paper II that even the Friedmann models do not seem to 
remain homogeneous along the past null geodesic (see figures 1, 2, 6 and 7), if the 
homogeneity of the standard model does not survive, where is the strong need for a 
crossover? In paper I it was assumed that the Tolman metric would eventually join 
a Friedmann background and, among other things, I argued the need for that was to 
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make the model compatible with a different interpretation of the Copernican principle. 
However, in the light of inhomogeneity even in the Friedmann metric, it could be argued 
that from an observational point of view, that is, in calculations along our past null cone, 
if the relativistic fractal cosmology developed in this series of papers has or has not a 
Friedmann background might well be irrelevant. 

7 Conclusion 

In this work I have presented numerical solutions of particular Tolman models of elliptic, 
parabolic and hyperbolic types featuring fractal behaviour along the past light cone. The 
initial conditions of the numerical problem is discussed in detail and I found three differ- 
ent sets of initial values that can, in principle, be used to solve the problem numerically. 
In practice, however, I have made use of only two due to strong stiff behaviour in the 
numerical solutions. The algorithm used to get the solutions has been described, as well 
as some details of the manner it was implemented in the computing code. 

The spatially homogeneous Friedmann models were treated as special cases of the 
Tolman solution and solved numerically. In an extension of the results of paper II, I 
have found that the K — ±1 Friedmann cases also do not appear to remain homoge- 
neous along the backward null cone, and this is explained because the models are spatially 
homogeneous: when the density is averaged along the null geodesic, going through hy- 
persurfaces of constant t, but with different values for the local density, this average 
changes. Furthermore, as in the open and closed Friedmann models the volume density 
decreases the farther we go, it would seem that similarly to the Einstein-de Sitter model 
studied in paper II, all cases of the standard cosmological model have zero global density 
at the big bang singularity hypersurface. 

Fractal solutions were obtained in all classes of Tolman's solution and the mathemat- 
ical criterion for getting them was their obedience to the de Vaucouleurs' density power 
law. Nonetheless, considering that observationally realistic models have to follow the 
linearity of the luminosity distance-redshift diagram, have a Hubble constant within the 
currently accepted uncertainty and fractal dimension in agreement with measurements 
of the 2-point correlation function, I found only hyperbolic type solutions obeying all 
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these criteria. 

The simplest fractal hyperbolic solution has simultaneous big bang and it was found 
that this model has fractal features in many different epochs, does not have homothetic 
self-similarity (although some of the parabolic fractal solutions do) and if we assume a 
Friedmann background spacetime, the matching between these two metrics imply that 
the background must be an open Friedmann model. A different scenario where the 
"Great Wall" could be modelled was also considered. The paper ends with a discussion 
on some objections raised by some authors against a fractal cosmology where I have 
shown how some are imprecise and inappropriate, leading to untenable objections to 
this sort of cosmology. 

And what about the origin of this fractal system? If a Friedmann background is 
assumed, it might be a cellular-type structure (Feng, Mo & Ruffini 1991; Fabbri & Ruffini 
1992 and references therein), but with the difference that an upper cutoff to homogeneity 
in the past light cone is not compulsory, as discussed in the previous section. In this 
case this structure could be a result of a fragmentation process (Ruffini, Song & Taraglio 
1988). If we do not assume a Friedmann background, Bonnor (1974) showed that a 
hyperbolic Tolman model can start inhomogeneous, remain inhomogeneous during the 
course of its evolution and at large times, when the dust has already escaped from its 
own gravitational field, there is nothing that could alter the model's density distribution. 
In this case the function f3(r) does not appear neither in the function R(r,t) nor in the 
density p, which was interpreted as meaning that at large t the model has "forgotten" 
its big bang, at least as far as the density is concerned. In conclusion, the origin of the 
fractal structure appears to be an open problem. 
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Figure Captions 

Figure 1: This graph is the same appearing in paper II and shows the analytical 
results of a log p v vs. log di plot for the Einstein-de Sitter model in the 
range 0.001 < r < 1.5 and with (3 = 2.7. The inhomogeneity of the 
model along the past null geodesic is clearly visible. 

Figure 2: Plot of the numerical results of the volume density p v against the luminosity 
distance di in the range < r < 1.5 for the Einstein-de Sitter model. Only 
steps bigger than 0.0025 are shown here. Distance is given in Gpc (10 9 
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Figure 3: 
Figure 4: 
Figure 5: 



Figure 6: 



Figure 7: 
Figure 8: 



pc) and mass in units of 2.09 x 1O 22 M . By comparing these results with 
the analytical ones shown in figure 1 we can see that the numerical scheme 
produces the expected results. 

Numerical integration of the null geodesic in the Einstein-de Sitter model. 
Time is in units of 3.26 Gyr. 



Numerical integration of the equation (19) which gives the function I(r) in 
the Einstein-de Sitter model. 

The same results of figure 2 for the Einstein-de Sitter model, but with all 
steps visible. The numerical code takes an excessive number of steps around 
the origin to integrate this model, despite the smoothness of the integrated 
functions t(r) and J(r) shown in figures 3 and 4, respectively. The solution 



in this graph shows clearly the stiff behaviour of equations (11) and (19) at 
around r = in this model. 

Numerical results of logp^ versus logrf/ in an open Friedmann model. Very 
small steps are omitted and so stiffness is not visible. From now on the 
figures will not have visible stiff behaviour to avoid saturation. 

Results of p v plotted against di for the integrations of a recollapsing Fried- 
mann model. 



Numerical results of the parabolic model obtained with functions (38) in 



tegrated in the interval < r < 2 and its straight line fitting carried out 
according to equation (37). The constants took the values a = 1CT 5 , p = 
0.9, Po = 2, t]q = 50, q = 0.65 and the resulting fitting coefficients are 
ai = —5.5 and 02 = —1.6. Considering equation (25) the fitting coefficients 
tell us that the fractal dimension of the dust in this model is D = 1.4 and 
the lower cutoff constant is o = 2.8 x 10 6 . 

Figure 9: Results of p v and di for the integration of the parabolic model with functions 
in the interval < r < 5. The constants used had the values: a = 
p = 0.5, p = 1.5, rji = 1000. The fitted straight line is also visible 
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and the resulting fitting coefficients are a± = —3.9 and ai = —1.3, which 
means a fractal dimension D = 1.7 and a lower cutoff constant a = 9.7x 10 7 . 

Figure 10: Distance-redshift diagram obtained by integrating the parabolic model of 



functions (38). Note the absence of a linear relation between the luminosity 
distance di and the redshift z. This fact makes this model observationally 
unrealistic. 

Figure 11: Distance-redshift diagram obtained by integrating the parabolic model of 



functions (39). The diagram is only partially linear, but even so the Hubble 



constant associated with this part is just at or below the lower limit for Hq. 



This point implies that functions (39) are not a good modelling for the 
observations. 

Figure 12: Volume density p v vs. luminosity distance di for the integration of the 



elliptic model given by functions (40) in the interval 0.001 < r < 0.07. The 



constants took the values a — 10, p = 1.4, j3 = 0.7, 770 = 50, q = 0.65 and 
the resulting fitting coefficients are a\ = —2.7 and 02 = —1.3, which means 
D = 1.7 for this model. 

Figure 13: Volume density p v vs. luminosity distance d\ for the integration of the 



elliptic model given by functions (41) in the interval 0.001 < r < 0.07. The 



constants used had the values: a = 10, p = 1.4, (3q = 0.7, r/i = 1000. The 
fitting coefficients obtained are a x = —2.4 and a<i = —1. The models has 
fractal dimension D = 2. 

Figure 14: Distance-redshift diagram obtained by integrating the elliptic fractal model 



produced by functions (40). The diagram is not linear which makes the 



model incompatible with observations. 
Figure 15: Distance-redshift diagram obtained by integrating the elliptic fractal model 



produced by functions (41). The diagram can be approximated to a linear 



relation, but that produces a too low value for H 
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Figure 16: 



Figure 18: 



Plot of the results of p v vs. di obtained by the numerical integration 
in the interval < r < 0.07 of the hyperbolic model constructed with 



functions (42). The values of the constants used to find these results are 
a = 10~ 4 , p = 1.9, Po = 3.6, 770 = 50, q = 0.65 and the fitting coefficients 
obtained are a\ = —7.3 and 02 = —1.6. 

Figure 17: Distance-redshift diagram obtained with the hyperbolic model of functions 



(42). The diagram is approximately linear but the associated Hubble con- 



stant of about 23 km/s/Mpc is too low comparing with the current uncer- 
tainty of Hq. 

Plot of the results of p v vs. di obtained by the numerical integration in the 
interval < r < 0.07 of the hyperbolic model constructed with functions 
(43). The constants used to find these results are a = 10~ 4 , p = 1.4, @ = 
3.6, 771 = 1000. The fitting coefficients obtained are a% = —6.2 and a 2 = 
— 1.7. These coefficients mean a fractal dimension D = 1.3 and a lower 
cutoff constant a = 5.4 x 10 5 . 

Figure 19: Distance-redshift diagram obtained with the hyperbolic model consisted of 



functions (43). This is one of the best linear approximation for a d\ vs. z 



diagram and the slope gives H = 61 km/s/Mpc, a value within the current 
uncertainty in the Hubble constant. 

Figure 20: Results of the volume density p v vs. luminosity distance di of the hyperbolic 



model (44) with simultaneous big bang. The integration is in the interval 
< r < 0.07 and the constants are a = 10~ 4 , p = 1.4, j3 = 3.6. The 
fitting coefficients calculated are a\ = —6.0 and a 2 = —1.6, giving a fractal 
dimension D = 1.4 and a lower cutoff constant a = 8.7 x 10 5 . 

Figure 21: Distance-redshift diagram obtained with the hyperbolic model of functions 



(44). This is the best diagram obtained for di vs. z and the slope gives 



Hq = 80 km/s/Mpc, a value which is not only within the current uncertainty 
in the Hubble constant but also agrees with recent measurements on it (see 
Peacock 1991). 
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Figure 22: Plot of the results for the cumulative number counting N c vs. the redshift 



z given by the integration of the model (44) with the same parameters as 
in figure 20. 
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